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Abstract 

We propose two algorithms for sampling from two gamma variates 
possessing a negative correlation. The case of positive correlation is 
easily solved, so we just mention it. The main problem is the lowest 
value of the correlation coefficient that can be reached. The starting 
point of both algorithms is generation from a bivariate density with 
uniform negatively correlated marginals. Actually the first method 
uses a degenerate bivariate density since it considers two uniforms re- 
lated by a linear relationship. Then we resort essentially to the inverse 
transform method. For both algorithms we stress restrictions on the 
parameters and rigidities. 

1 Introduction 

To fix notation, if X follows a gamma distribution with parameters A and 
a, that is X ~ ^(A, a), then the density is 

^Og— Ax^o— 1 

fx{x) = z— , X > 0. 

r{a) 

Johnson and Kotz 1972] define a multivariate gamma distribution in the 
following way. We have m + 1 independent standard (that is, each with 
A = 1), gamma variates (in general with different a's) Xq, Xi, Xm- 
Define 

Yj = Xo + Xj, j = 1, m 
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Then the Yj, j = 1, . . . , m are distributed according a m-variate gamma 
variable. We can see that each marginal Yj is a standard gamma variable 
with a = ao + aj and 



Cov{Yj, Yj, 



Var(Xo) 



and so all the marginals are positively correlated. 

We follow this suggestion, with m = 2, but we exploit the preserving mono- 
tonicity property of the inverse transform method, see Fishman 1996] to 
generate a couple of gamma variates, X\ and X2, with a rigid amount of 
negative correlation and then generate Xq to allow for some flexibility. 

2 First Method 

Let X\ ~ ^(1, r) and X2 ~ ^(1, s), with r and s integers. Assume, without 
loss of generality, r < s. Let Ui be s independent uniforms on the unit 
interval, that is Ui ~ ^(0, 1), i = 1, . . . , s. Then we can write, using 
properties of the gamma variate and the inverse transform method. 



Xi = -5]lnC/i, 

i=\ 

s 

X2 = -Y.Hi-Ui). 



i=l 



It follows that 

Cr s 
^InC/i, ^ln(l - Ui 
i=l 1=1 

r 

= 5]Co?;(ln[/i, ln(l-C/i)) 




Now define, with Xq ~ G{1, ao) independent of Xi and X2, 

Yi=Xq + X^, 
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Y2=Xq + X2. 

It follows Yi ~ ^(1, ckq + ^) aiid I2 ~ ^(1, + s). Furthermore 
Cot;(yi, ^2) = Var(Xo) + Co?;(Xi, X2) 
= ao + 
Because 

c = 1 = -0.644934, 

6 

we can see that varying r and ao we have some freedom in modelling a 
negative covariance. 

The correlation coefficient between Yi and I2 is given by 

p{Yi, Y2) = ==. (1) 

V(ao + r){ao + s) 

For oq + rc < 0, which guarantees a negative correlation, we see that 

dp{Yi, Y2) 



>0, 



and so when the correlation is negative, it is an increasing function of Oq. 
The lower bound for p is given by 



rc ir 



y/rs V s 

When r = s and assuming that oq = means that Xq = with probability 1 
we reach the most negative correlation possible between two gamma variates. 
If we are given + r = m, + s = n and p = po then from Equation |l] we 
get 

m + r(c- 1) 

Po = ^= , (2) 

^/nln 



which can be solved for r and subsequently obtaining oq and s. The draw- 
back is that r and s are required to be integers, so we cannot arbitrarily 
choose m, n and p. To give an idea of the situation we present some tables, 
assuming for simplicity that both m and n are integers. It follows that also 
ao has to be an integer. Because we want a negative correlation and because 
ao > we have the restriction 

r < m < r(l — c). 
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Using Equation ^ we obtain 



r 


TO 


s 


P 


r 


TO 


s 


P 


r 


m 


s 


P 


2 


2 


3 


-.5266 


2 


2 


5 


-.4078 


2 


2 


8 


-.3224 


2 


3 


4 


-.0837 


2 


3 


6 


-.0683 


2 


3 


9 


-.0557 


5 


5 


6 


-.5887 


5 


5 


8 


-.5098 


5 


5 


11 


-.4348 


5 


6 


7 


-.3432 


5 


6 


9 


-.3027 


5 


6 


12 


-.2621 


5 


7 


8 


-.1636 


5 


7 


10 


-.1463 


5 


7 


13 


-.1283 


5 


8 


9 


-.0264 


5 


8 


11 


-.0239 


5 


8 


14 


-.0212 


8 


8 


9 


-.6080 


8 


8 


11 


-.5500 


8 


8 


14 


-.4875 


8 


9 


10 


-.4384 


8 


9 


12 


-.4002 


8 


9 


15 


-.3579 


8 


10 


11 


-.3012 


8 


10 


13 


-.2771 


8 


10 


16 


-.2497 


8 


11 


12 


-.1879 


8 


11 


14 


-.1740 


8 


11 


17 


-.1579 


8 


12 


13 


-.0928 


8 


12 


15 


-.0864 


8 


12 


18 


-.0788 


8 


13 


14 


-.0118 


8 


13 


16 


-.0110 


8 


13 


19 


-.0101 


12 


12 


13 


-.6196 


12 


12 


15 


-.5768 


12 


12 


18 


-.5265 


12 


13 


14 


-.4995 


12 


13 


16 


-.4672 


12 


13 


19 


-.4288 


12 


14 


15 


-.3960 


12 


14 


17 


-.3720 


12 


14 


20 


-.3429 


12 


15 


16 


-.3059 


12 


15 


18 


-.2884 


12 


15 


21 


-.2670 


12 


16 


17 


-.2267 


12 


16 


19 


-.2144 


12 


16 


22 


-.1993 


12 


17 


18 


-.1565 


12 


17 


20 


-.1485 


12 


17 


23 


-.1385 


12 


18 


19 


-.0940 


12 


18 


21 


-.0894 


12 


18 


24 


-.0836 


12 


19 


20 


-.0379 


12 


19 


22 


-.0361 


12 


19 


25 


-.0338 



Nothing essentially changes if we require gamma's with the first parameter 
different from 1: it is enough to multiply Yi and I2 by, say, a constant /?. 
The correlation coefficient is of course not affected. 

3 Joint Density 

Obtaining the joint density of Yi and Y2 in the general case is cumbersome. 
To give an idea we consider the simplest case, r = s = 1, and we set a < —c 
that is Q < 0.644934, to guarantee a negative correlation. So we have 

Yi = Xo-lnU, (3) 

Y2 = Xo-Hl-U). (4) 
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with independence of Xq and U. It turns out that marginahy Yi and I2 are 
both g{l, l + ao). Setting 

[2/2 = xq- ln(l - u), 

we get 

C 1 



\ XQ = yi - In (1 + 62^1 -2^2] 
Because > we have the inequahty 

yi - In (1 + e^^-y^) > 0. 

The solution to the equation 

yi-ln{l + ey^-y^) =0 

is given by 

2/2 = yi - In (e^i - 1) . 
It follows that we have xq > if 

y2>yi-ln(e2'i-l). 

For the Jacobian of this transformation we have 

p2/l-J/2 

J = 



+ 

The joint density of Xq and U is given by 

fxo,u{xo, u) =e-^°x^"\ 
The quantity e~^° J simplifies to 



< n < 1 

< Xq. 



and consequently we get for the joint density of Yi and Y2 
[2/i-ln(l + efi-2'2)pi f 2/1 > 



evi +ey^ ' I y2 > yi - in (e^i - 1) . 
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We see that for a = 1 this simphfies to 

fnMvu y2) = -iJT^,, I > -ln(e^^ - 1) 

but in this case the correlation is positive: 

1 + c 



2 

0.177533. 



In this last situation the marginals are ^(1, 2). This can be checked using 
the fact that 

1 -6 + In fa + e'') 
■ dx = . 



a + e^ 

4 Second Method 

Again looking for flexibility we analyze another method where we start gen- 
erating s samples (C/ij, from a bivariate distribution with density 

/(ui, ^2) = 1 + 0(1 -2ni)(l-2n2), < m, ^2 < 1. (5) 

This density is of the form studied by Long and Krzysztofowicw 1995]. 
For — 1 < < 1 this function is a proper density. It turns out that marginally 
C/i- and {72i ai^e uniforms over the unit interval and the correlation coefficient 
is given by 



Now define 



3 

r 

i=l 
s 

X2 = -^lnC/2,. 



2=1 



Now 



Cou(lnC/i,, ln;72j = f\nxlny{l + 0{l - 2x){l - 2y)) dxdy - I 

Jo Jo 
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It follows that 



CoviXi, X2 



rO 



Proceeding as in the other method we define, with Xq ~ t/(I, uq) indepen- 
dent of Xi and X2 , 



Y2 = Xo + X2 



and we get 



"0 + — 



(6) 



V(ao + r){aQ + s) 

Now we have a negative correlation if AaQ+rO < 0. Because of this inequality 
the admissible range for 9 is —1 < 6 < 0. The lower bound for p is 

e [¥ 

so we cannot hope to do better than 

p > -0.25. 

Repeating the same reasoning that led us to Equation ^ now we have 

4m - r(4 - 9) 



Po 



4\/mn 



The restrictions are now 



Solving for r we get 



r < m < 



r(4 - 9) 



4m — ApQ^mn 



(7) 



(8) 



Other restrictions on the minimum value admissible for po arise from Equa- 
tion ^ and the fact that r must be integer. The minimum value of r as given 
by Equation is obtained for 9 = —1. Then 



4m — ApQ^/mn 



< m — 1, 



which implies 



Po > 



m — b 
AJmn 



(9) 



7 



Under the condition m > 6 this attainable lower bound is an increasing func- 
tion of n and because n > m it is a decreasing function of m: in particular 
if m = n the limit for m going to infinity is —0.25. 

Once observed the conditions on m and p we have the following algorithm. 



1. Set y = Am — ApQyJran. 



2. Obtain a = |. 
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3. Set r = [a] , that is the lowest integer not lower than a. Because of 
the construction such an r exists and r < m. 

4. Obtain 6 = A - y/r. 

As an example, imagine m = 7, n = 10. The attainable lower bound is 
-0.0597. Set for example po = -0.05. Then | = 5.2653. Take r = 6, so 
that 6* = 4 - 4.38778 = -0.38778. 

Once is obtained, the next step is to evaluate oq = m — r and s = n — oq. 
Then we generate s samples from the density given in Equation ^: this can 
be done for instance with the acceptance-rejection method. And then we 
follow the final construction to obtain Yi and Y2. 

We can note that in this second method the attainable lower bound for p is 
sensibly greater than in the first method, but for the admissible values we 
have a complete flexibility. The other negative point is that we have to gen- 
erate samples from a bivariate density, whose covariance scaler 9 (as termed 
by Long and Krzysztofowicw) has to be preliminarily obtained, instead of 
generations from univariate densities. 
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